Astronomy Sz Astrophysics manuscript no. ms3717 


February 2, 2008 


(DOI: will be inserted by hand later) 





The mass-velocity and intensity-velocity relations in jet-driven 

molecular outflows 

Turlough P. Dowries 1 and Sylvie Cabrit 2 



O 
O 



> 

ON 

(N 
m 
m 
O 
m 
o 

Oh 
6 

C3 



13 



School of Mathematical Sciences, Dublin City University, Dublin 9, Ireland 
2 LERMA, Observatoire de Paris, 61 Av. de l'Observatoire, F-75014 Paris 

Received June 18, 2001;accepted March 12, 2003 

Abstract. We use numerical simulations to examine the mass-velocity and intensity-velocity relations in the CO 
J=2-l and H2 S(l)l-0 lines for jet-driven molecular outflows. Contrary to previous expectations, we find that the 
mass- velocity relation for the swept-up gas is a single power-law, with a shallow slope ~ —1.5 and no break to a 
steeper slope at high velocities. An analytic bowshock model with no post-shock mixing is shown to reproduce 
this behaviour very well. 

We show that molecular dissociation and the temperature dependence of the line emissivity are both critical in 
defining the shape of the line profiles at velocities above ~ 20 km s . In particular, the simulated CO J=2-l 
intensity- velocity relation does show a break in slope, even though the underlying mass distribution does not. 
These predicted CO profiles are found to compare remarkably well with observations of molecular outflows, both 
in terms of the slopes at low and high velocities and in terms of the range of break velocities at which the change 
in slope occurs. Shallower slopes are predicted at high velocity in higher excitation lines, such as H2 S(l)l-0. 
This work indicates that, in jet-driven outflows, the CO J=2-l intensity profile reflects the slope of the underlying 
mass- velocity distribution only at velocities < 20 kms , and that higher temperature tracers are required to 
probe the mass distribution at higher speed. 
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1. Introduction 

' Various authors have noted that the intensity-velocity re- 
lationship observed in low-J CO lines in molecular out- 
flows tends to follow a broken power-law Ico{v) 0(1 u_7 j 
with 7 w 1.8 ± 0.5 up to line-of-sight velocities Vbrcak ~ 
' 10-30 kms -1 and 7 « 3-7 at higher velocities (e.g. 
Rodriguez et al. Stahler Bachiller & Tafalla 

1999] Richer et al. 2000J. This property is an important 
test for proposed mechanisms of molecular outflow accel- 
eration. In particular, recent work has addressed this issue 
in the case of entrainment by a jet. 

Using an analytic model of a jet/bowshock system, 



Zhang & Zheng (1997) predicted a mass distribution, 
m(v), following a broken power-law with slopes u. « 
1.8 up to 10 kms -1 and /i w 5.6 beyond. The observed val- 
ues of 7 in CO intensity profiles are then reproduced, pro- 
vided the CO abundance and line excitation do not vary 
much with velocity. However, this seems highly unlikely, 
given the broad range of shock strengths in a bowshock. 
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Smith, Suttner & Yorke i 1 997 i conducted numerical 
simulations of jet-driven molecular outflows that take into 
account dissociation and heating in shocks. They find that 
Ico{ v ) follows a power-law i; -7 with 7 ssl.2— 1.6 up to 
10 kms -1 , and a steeper slope further out, in qualitative 
agreement with observations. They attribute this result 
to a much steeper mass- velocity relation, m(v) oc u~ 3 - 5 , 
than Zhang & Zheng l|1997fl . However, their reasoning in- 
volves an erroneous high-temperature dependence of the 
CO emissivity (oc T instead of T -1 ). Therefore, the ac- 
tual origin of the slope of Ico{v) in jet simulations, and 
the underlying mass-velocity relation itself, remain to be 
clearly established. In this work we clarify this issue using 
simulations at higher resolution, and analytical modelling. 

2. Numerical method 

The code used in this work is very similar to that used 
in Downes & Ray (1999). The initial conditions are also 
very similar, but we give a brief overview of them here for 
completeness. 

The simulations are performed in 2D cylindrical sym- 
metry. The densities of molecular and atomic hydrogen 
are tracked, along with the ionisation fraction of hydrogen. 
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The CO density is assumed to be a constant (10 -4 ) frac- 
tion of the H2 density by number. The numerical scheme 
is a Godunov scheme which is second order in time and 
space (see, e.g., Downes & Rav ll999"|l . 

The grid-spacing was set at 10 14 cm. The jet radius 
was 5 x 10 15 cm, and its density was 100 cm -3 (equal to 
the ambient density), thus allowing a resolution of cool- 
ing layers a factor of 10 2 -10 4 higher than in Smith et al. 
(1997). The time-averaged jet velocity, vq, was 215 km 
s . Superimposed on this were sinusoidal variations with 
periods of 5, 10, 20 and 50 years with total amplitude v\. 

In order to explore the effects of varying the jet veloc- 
ity, three simulations were run, as follows: (1) a "steady" 
jet where ^ = 0; (2) a "pulsed" jet where ^ = 0.6; and 
(3) a "ramped-up" jet where ^ = 0.6 and vq increases 
linearly in time from to 215 km s -1 over 100 years. 
Typically, simulations of jets start off with the jet propa- 
gating at full speed into the ambient medium. It is likely 
that YSO jets do not 'switch on' impulsively, and it has 
been suggested recently (Lim et al. l2002f) that a slow start- 
up of a jet may lead to increased molecular abundance at 
the head of the bowshock, which would then affect the re- 
sulting intensity-velocity relation that we wish to model. 
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Fig. 2. Plots of the mass- velocity relations for all swept- 
up material (i.e. excluding jet material), m(v), and for 
molecular swept-up material only, mH 2 (ti), for the pulsed 
jet at t = 400 yrs. Also shown are the intensity-velocity 
relations for the CO J=2-l line and for the H2 S(l)l-0 line 
(arbitrary vertical offsets). An angle of 30° to the plane of 
the sky is assumed. 



3. Numerical results 

In this section, we concentrate on the results of simulation 
(2) (pulsed jet) described in Sect. [3 a t a time t = 400 yrs. 
Results for other times and for the steady and ramped-up 
jets (simulations (1) and (3) in Sect.[2J) will be discussed 
in Sect.JSJ 

Figure ^ contains a log-scale plot of the density dis- 
tribution for the chosen simulation. We can see that the 
bowshock is rather irregular. These irregularities are at- 
tributed to the growth of the Vishniac instability (Downes 
& Ray 1999). In addition, several "mini-bowshocks" are 
present along the jet length, tracing internal working 
surfaces resulting from the variability of the jet veloc- 
ity. Figure [21 contains plots of various mass- velocity and 
intensity-velocity relations for the same simulation, as- 
suming an inclination of 30° to the plane of the sky. We 
discuss each of these distributions in detail in the follow- 
ing. 

3.1. The total mass-velocity relation 

A significant feature of Fig. [21 is that m(v), the mass- 
velocity relation for all swept-up material (i.e. excluding 
jet material), remains quite shallow across the whole ve- 
locity range. It is essentially flat at the lowest velocities 
(v < 3 kms -1 ), then follows an approximate power-law 
m(v) oc t> -M with exponent /1 ~ 1.5 at intermediate ve- 
locities, before rising slightly again above v ~ 40 kms -1 . 
The power-law slope that we find agrees quite well with 
the value of /1 ~ 1.8 predicted by the analytical model of 
Zhang & Zheng (1997) for low velocities, but their pre- 
dicted break to a steeper slope n ~ 5.6 at velocities above 
10 kms -1 is not seen. 



3.2. The molecular mass-velocity relation 

The mass-velocity relation for swept-up molecular mate- 
rial, mHj(«), behaves similarly to the total mass- velocity 
relation, m(v), at low velocities. It does, however, become 
considerably steeper than m(v) at higher velocities (v > 
20 kms -1 ). This steepening occurs because material at 
these velocities has been accelerated by faster shocks near 
the apex of the bowshock, where significant molecular dis- 
sociation occurs. Indeed, the typical dissociation limit for 
low-density non-magnetic shocks is ~ 30 kms -1 (see e.g. 
Flower et al. 2003). We note that the molecular fraction 
appears to rise again at the highest velocities (above 80 
kms -1 ). This is an unavoidable result of numerical diffu- 
sion at the bow head, which results in mixing of molecular 
jet material into the high-velocity swept-up gas. 

3.3. The Ico(v) relation 

The CO J=2-l line emissivity at each grid point of the 
simulation was calculated from the local density and tem- 
perature according to the analytical formulae of McKee et 
al. I|1982|l . which take into account sub-thermalization of 
the levels. Emission was assumed optically thin, which is 
probably a good approximation for velocities above a few 
kms -1 . 

Figure [21 shows that the resulting intensity-velocity re- 
lation for the CO 3—2-1 line, Ico( v )i follows the same 
slope as toh 2 (v) at low velocities, but breaks to an even 
steeper relation than toh 2 {v) at velocities above ~ 30 
kms -1 . This steepening results from the temperature de- 
pendence of the CO(2-l) emissivity per molecule. This 
dependence (assuming LTE) is given in terms of a func- 




tion fl(T) in Cabrit & Bertout (1990), and is plotted as a 
function of temperature for the CO J=2-l line in Fig. 13 
of Lada & Fich Ijl99()|l . The line emissivity per molecule is 
seen to initially increase steeply with T up to a maximum 
at T ~ 20 K, similar to the excitation energy of the upper 
level of the transition. For higher temperatures, however, 
it then decreases as T _1 , due to the larger number of en- 
ergy levels available to the molecule (or equivalently, as an 
effect of the partition function). Hence, as we go to higher 
velocities, not only are there fewer molecules to emit (due 
to shock dissociation), but they are hotter (because they 
have been through a stronger shock) and so they emit less 
efficiently in the CO J=2-l line. This produces an even 
steeper slope at high velocities in ico(^) than in mn 2 (v). 

3.4. The Ir 2 (v) relation 

It is interesting to explore the predicted intensity-velocity 
relation for the H2 S(l) 1-0 line, In 2 (v), as it has been 
recently observed in a few flows (Salas & Cruz-Gonzalez 
2002) and the line is of much higher excitation than the 
CO J=2-l line. Such predictions from jet simulations have, 
to the best of our knowledge, never been presented so far. 
For simplicity, we computed the line emissivity assuming 
LTE. Although this is not fully justified at the moderate 
densities of our simulations, it is sufficient to illustrate the 
change of slope with excitation of the line that we predict 
to occur in jet-driven outflows. 

From Fig. [5] we see that, unlike the CO 2-1 line, this 
intensity- velocity relation is shallower than mn 2 (v) up to 
30 kms -1 , and therefore does not show as dramatic a 
break in slope at high velocities. This is due to the fact 
that the the upper energy level for the H2 S(l) 1-0 line 
lies at 7000 K; hence the emissivity per molecule is still 
increasing with T up to temperatures of ~ 10 4 K, so that 
higher velocity, hotter material now emits more efficiently 
in the H2 S(l) 1-0 line. The In 2 (v) distribution is then 
shallower than the underlying mass distribution. 



4. Comparison with observations 

In this section, we compare our simulation results directly 
with observations. Figure 01 contains plots of our results 
for the m(v) and Ico{v) relations, along with observed 




log(v) (km/s) 

Fig. 3. Plots of the observed intensity-velocity relations 
for L1448, Orion A, NGC2071, L1551 and Mon R2 ('+' 
signs, crosses, stars, open boxes, and filled boxes respec- 
tively). Also shown are the simulated m(v) and Ico(v) 
(dotted and solid lines, respectively, with arbitrary ver- 
tical offsets) for comparison. The fits to both L1448 and 
Orion A assume an angle of 60° to the plane of the sky, for 
NGC2071 the angle is assumed to be 30°, and for L1551 
and Mon R2 the flows are assumed to be in the plane of 
the sky. 



intensity-velocity relations in CO J=2-l for various bipo- 
lar outflows (taken from Bachiller & Tafalla 1999). A sim- 
ilar comparison with the results of Salas & Cruz-Gonzalez 
( 2002) will be made in a forthcoming paper, using NLTE 
calculations of the H2 S(l) 1-0 line intensity. 

It is notable that the simulation results for the Ico (v) 
relation reproduce the observations remarkably well, given 
that the velocity of the jet was arbitrarily chosen and the 
simulation timescale is shorter than that of the physical 
systems (estimated ages range from 1000 yrs in Orion to 
several times 10 4 yrs for L 1551 and Mon R2). It is also 
interesting to note that the predicted Ico (v) more closely 
matches the observations than the total mass- velocity re- 
lation m(v), save for L1448. 

The only parameters which have been tuned in this 
comparison are the vertical scale (i.e. a scaling in the ambi- 
ent gas density) , and the assumed inclination to the plane 
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Fig. 4. Plots of the mass- velocity distribution of swept-up 
material for three simulations with differing jet velocity 
variability at t = 400 yrs, assuming an angle of 30° to 
the plane of the sky. Note the similarity between all three 
relations. 



of the sky. The latter parameter can be seen to have a 
sizeable effect on the break velocity. While our adopted 
inclinations do not exactly match that inferred from ob- 
servations for each of the systems shown, the trend is cor- 
rect in the sense that Orion A and L1448 are the most 
inclined to the plane of the sky, NGC2071 is intermediate, 
and L1551 and Mon R2 are believed to be closest to the 
plane of the sky. 

In summary, it seems clear, on the basis of these simu- 
lations, that a jet-driven bowshock can reproduce the ob- 
served Ico (v) remarkably well, both in terms of the slopes 
at low and high velocities and in terms of the range of 
break velocities where this change of slope occurs, without 
the need to invoke two distinct entrainment mechanisms 
as is done in Zhang & Zheng (1997). 

The simulation results do tend to be slightly too flat at 
the lowest velocities, especially in the oldest flows (Mon R2 
in particular). This could be due to the short timescales 
in our simulations, as we now show. 

5. The effect of age and time-variability 

5.1. Effect of time-variability of the jet velocity 

It is important to verify that our results are not criti- 
cally dependent upon the adopted jet velocity behaviour. 
Figure^Jcompares the m(v) relation for the pulsed jet with 
the swept-up mass- velocity relation for the other two sim- 
ulations (steady jet and ramped-up jet), at the same time 
and same inclination. It is clear that the relations are very 
similar. This is also the case for the mn 2 (v), Ico(v) and 
In 2 (v) relations (not shown). This strongly suggests that 
variability of the jet velocity (at least on time-scales less 
than about 50 yrs) will not greatly affect the resulting 
mass- velocity or intensity- velocity relations. 
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Fig. 5. Plots of the mass- velocity relation for all swept-up 
material in the pulsed jet (simulation (2)) at t = 100, 200, 
300 and 400 yrs, assuming an angle of 30° to the plane of 
the sky. 

5.2. Effect of age 

It is also important to verify that our results do not de- 
pend critically on the simulation timescale, since observed 
flows cover a large range of ages. Figure El contains plots 
of the mass- velocity relation for the pulsed jet simulation 
(2) at different times from t = 100-400 yrs. It can be 
seen that the part of the m(y) distribution that follows 
a power-law oc i>~ 1,5 extends to progressively lower ve- 
locities as time goes, while the flat portion of m(v) at 
very low velocity shrinks. This behaviour can be under- 
stood as follows: the power-law regime corresponds to the 
part of the bow structure presently on the grid, while the 
flat part of m(v) comes from the truncation of the bow 
at the left border of the computational grid; as the bow 
head advances, more of the bow structure appears on the 
grid, so that the part affected by truncation is at a lower 
speed (since velocity decreases away from the bow head). 
Hence, we expect that for longer timescales, the flat part 
of m(v) and Ico(v) w ih fell to very low velocities, while 
the ii -1 - 5 power-law should extend over the whole observ- 
able velocity range, in better agreement with observations 
of old flows such as Mon R2. Simulations over very long 
timescales are under way to study this hypothesis. 

6. Analytic model 

In this section we explore in more depth how a jet-driven 
bowshock can produce the power-law slope fj, ~ 1.5 of 
the mass- velocity relation discussed in Sect. |3~T1 Figure [B] 
contains a schematic diagram of the set-up used in this 
discussion. 

Consider an idealised bowshock, the shape of which is 
given by the relation 2 oc r s , moving in the positive z di- 
rection with constant velocity vq into an ambient medium 
of constant density p a . We assume that there is no mixing 
in the bow, so that the magnitude of the post-shock veloc- 
ity in the observer's frame is vqSuiO, with tan 9 — dr/dz 
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Fig. 6. Schematic diagram of the set-up used in the ana- 
lytic calculations of the mass- velocity relation. 

(see Fig. EJ. An observer situated along the z-axis will 
then see a line-of-sight velocity of 



Vq sin 



(1) 



Far from the apex, where 9 is small, v oc tan 2 9 oc (r 2 ) 1 s . 
The rate of increase of mass at projected velocity v is then 

— (m(v)dv) = p a vo2Trr(v)dr oc (u)^ 1 s ^ dv. (2) 

Hence the mass-velocity relation in this case is a power- 
law 



m(v) oc v M with /i 



s - 1 



(3) 



In order to verify that the no-mixing model can in- 
deed explain the m(v) relations found in our simulations, 
we need to determine the appropriate value of s for the 
simulated bows. For this, we examine the rate at which 
ambient mass is swept up by the bowshock, integrated 
over all velocities. For an adopted ideal bow shape z oc r s , 
this rate is 



dm 
~dt 



(t) = wr%Pa.v 



(4) 



where tq oc z l ' s (t) is the maximum bowshock radius and 
z(t) = vot is the length of the bowshock, both at time t. 
Hence, the total swept-up mass as a function of time is 
predicted to vary as 



m(t)oct a with a = l + 2/s. 



(5) 



In Fig. \7\ we plot the total swept-up mass m(t) as a 
function of time for each of the simulations. Indeed, we 
find that it closely follows a power-law, with a slope a = 
1.7-1.8. Note that the apparent deviation of simulation (3) 
from this law at early times is due to the initial ramping 
up of the jet velocity vq. Since the bow model predicts 
a = 1+2/s, we infer that the appropriate value of s for our 
simulated bowshocks is s = 2.5-2.9. If we now substitute 
this range of values of s into equation [21 we find that our 
simple bow model with no mixing predicts an exponent of 
the mass-velocity relation in the range /i = 1.5—1.7. This 
is remarkably close to what we find in the simulations 
(Sect. I3~TJ> . Note that the same bow model, but assuming 
instantaneous mixing of post-shock swept-up gas, predicts 
a much steeper slope fi = 2+s/2 = 3.5 (Smith et al. H997|l . 
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Fig. 7. Plots of the total swept-up mass against time for 
each of the simulations. 




Fig. 8. Plot of the total swept-up mass as a function of 
velocity predicted for an ideal bowshock shape described 
by z oc r 3 viewed pole-on, with no post-shock mixing. The 
straight line is the power law of exponent -1.5 predicted 
for the bow wings (see text). 



It is important to keep in mind that the power-law 
slope is predicted to hold only for small values of 9, i.e. 
sufficiently far from the head of the bowshock. In Figure 
El we plot the mass- velocity relation for a theoretical bow 
with s = 3, calculated numerically without the assumption 
of small 9. The power-law of slope /i = s/(s — 1) = 1.5 is 
seen to give a very good approximation to the actual m(v) 
up to about 20% of the jet velocity. At higher velocities, 
the mass-velocity relation flattens out and then increases 
with velocity. This behaviour is clearly also seen in our 
simulated m(v) (see Fig. [3J|, although we do not find as 
sharp a peak at vo, due to the fact that the head of the 
bowshock is rather irregular. 
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7. Conclusions 

We find that the mass-velocity relation in simulations of 
jet-driven molecular outflows is a shallow power-law of 
exponent ~ -1.5, with no break to a steeper slope at higher 
velocity, unlike previous analytical predictions by Zhang 
& Zheng ( 1997$. This m(v) relation can be well explained 
by a bowshock model with no post-shock mixing. It does 
not appear to be consistent with the full mixing model 
proposed by Smith et al. Q1997[) . 

We also find that the resulting intensity-velocity rela- 
tion for the CO J=2-l line compares remarkably well with 
observations of molecular outflows. In particular, it does 
have a break in slope around 20-30 km s _1 . The break is 
found to result from molecular dissociation near the bow 
apex, and from the 1/T-dependence of emission at tem- 
peratures exceeding the energy of the upper level of the 
line. Because of this dependence on T, a jet-driven model 
predicts a shallower slope at high velocity in higher exci- 
tation lines (e.g. H2 S(l)(l-0) and high-J CO lines), which 
could be tested by ongoing studies. Another implication 
of our results is that, in jet-driven outflows, the CO J=2-l 
intensity profile reflects the slope of the underlying mass- 
velocity distribution only at velocities < 20 kms -1 , and 
that higher temperature tracers are required to probe the 
mass distribution at higher speed. 

The only weak point of a bowshock model with no 
mixing is that emitting gas expands perpendicular to the 
bow surface, and thus there is some difficulty in explain- 
ing the apparent "forward-directed" motion of molecular 
outflows noted by Lada & Fich lfHJ 96l . Quantitatively, our 
simulations predict comparable amounts of redshifted and 
blueshifted gas up to velocities of 3-4 kms -1 . Full mix- 
ing was invoked by Smith et al. (1997) to alleviate this 
problem, but we have shown here that it would predict 
a very steep m(v) of slope ~ —3.5 and an even steeper 
CO intensity-velocity relation (due to dissociation and 
temperature effects), inconsistent with observations. We 
note that another situation in which motion will be more 
forward-directed is when the bowshock propagates into al- 
ready moving material. This is an interesting issue which 
clearly requires further work. 
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